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Abstract 



We report quantum Monte Carlo calculations of ground and low-lying ex- 
cited states for A = 8 nuclei using a realistic Hamiltonian containing the 
Argonne fig two-nucleon and Urbana IX three- nucleon potentials. The cal- 
culations begin with correlated eight-body wave functions that have a filled 
a-like core and four p-shell nucleons L5-coupled to the appropriate {J^;T) 
quantum numbers for the state of interest. After optimization, these varia- 
tional wave functions are used as input to a Green's function Monte Carlo 
calculation made with a new constrained path algorithm. We find that the 
Hamiltonian produces a '^Be ground state that is within 2 MeV of the exper- 
imental resonance, but the other eight-body energies are progressively worse 
as the neutron-proton asymmetry increases. The ^Li ground state is sta- 
ble against breakup into subclusters, but the ^He ground state is not. The 
excited state spectra are in fair agreement with experiment, with both the 
single-particle behavior of ^He and ^Li and the collective rotational behavior 
of ®Be being reproduced. We also examine energy differences in the T = 1 
and 2 isomultiplets and isospin-mixing matrix elements in the excited states 
of '^Be. Finally, we present densities, momentum distributions, and studies of 
the intrinsic shapes of these nuclei, with ^Be exhibiting a definite 2a cluster 
structure. 

PACS numbers: 21.10.-k, 21.45. -Fv, 21.60.Ka, 27.20. -Fn 
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I. INTRODUCTION 



The A = 8 nuclei have many interesting properties that we would like to understand 
on the basis of the bare interactions between individual nucleons. The three strong-stable 
nuclei, ®He, *Li, and ^B, decay by weak processes with half-lives of ^ 1 second down to 
^Be, which then immediately fissions to two ^He nuclei; thus there are no long-lived A = 8 
nuclei. These facts have tremendous consequences for the nature of the universe we live 
in: they make the production of elements beyond A = 7 in the early universe very difficult 
and help give stars like the sun a long stable lifetime. The ground and excited states of 
the A = 8 nuclei also provide a very sensitive testing ground for models of nuclear forces. 
In particular, ^He is the most neutron-rich strong-stable nucleus, and thus an ideal place 
to study the isospin dependence of the three- nucleon force. The T = 1, 2 isomultiplets are 
also a good place to look at charge-independence breaking, while some excited states in ^Be 
display significant isospin mixing. 

Previously we have reported variational Monte Carlo (VMC) and Green's function Monte 
Carlo (GFMC) calculations for the A = Q,7 nuclei In this paper we extend these cal- 
culations to the ground and excited states of A = 8 nuclei. (Some preliminary A = 8 
results were given in a number of conference proceedings We first construct a trial 

function, "^t, with the proper {J'^]T) quantum numbers and antisymmetry, and optimize 
the energy expectation value. We then use as the starting point for a GFMC cal- 
culation, which projects out the exact lowest-energy state by the Euclidean propagation 
\E'o = lim^^oo exp[— (if — Eq)t]'^t- We believe the resulting energy estimates are accurate 
to 1 to 2% of the binding energy in most cases. 

In our previous work, the Argonne fig plus Urbana IX Hamiltonian was fairly successful 
in generating the excitation spectra of the A = 6,7 nuclei, but did not give quite enough 
binding in the lithium isotopes, while ^He was unstable against breakup. One of our main 
goals in studying the A = 8 nuclei is to continue to test the Hamiltonian, and lay the ground 
work for studies of more sophisticated and accurate three-nucleon forces. The Hamiltonian is 
reviewed in Sec. II. The variational wave functions and calculations are described in Sec. HI. 
Our GFMC calculations make use of a new, more efficient algorithm called constrained-path 
propagation, which is described in Sec. IV. Results for the ground-state energies, excitation 
spectra, isomultiplet differences, and isospin-mixing matrix elements are given in Sec. V. 
In Sec. VI we present various density distributions, and in Sec. VII we discuss the intrinsic 
shapes of these nuclei. Our conclusions are given in Sec. VIII. 

II. HAMILTONIAN 

These calculations all use the same realistic Hamiltonian, which includes a nonrelativistic 
one-body kinetic energy, the Argonne f is two-nucleon potential and the Urbana IX three- 
nucleon potential @]: 

^ = E^^ + E^^.+ E (2-1) 

i i<j i<j<k 

The kinetic energy operator has charge-independent (CI) and charge-symmetry-breaking 
(CSB) components, the latter due to the difference in proton and neutron masses. 
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= Kf^ + Kf'^ ^ -^(— + — )V,^ - ^(— - —)r^Nl . (2.2) 

The Argonne fig model is one of a class of new, highly accurate NN potentials that fit 
both pp and np scattering data up to 350 MeV with a x^/datum near 1. The potential can 
be written as a sum of electromagnetic and one-pion-exchange terms and a shorter-range 
phenomenological part: 

% = ^^+<^+4- (2-3) 

The electromagnetic terms include one- and two-photon-exchange Coulomb interactions, 
vacuum polarization, Darwin-Foldy, and magnetic moment terms, with appropriate proton 
and neutron form factors. The one-pion-exchange part of the potential includes the charge- 
dependent (CD) terms due to the difference in neutral and charged pion masses. The shorter- 
range part has about 40 parameters which are adjusted to fit the pp and np scattering data, 
the deuteron binding energy, and also the nn scattering length. The one-pion-exchange and 
the remaining phenomenological part of the potential can be written as a sum of eighteen 
operators, 

+<= E ^pin,)Ot, . (2.4) 

p=l,18 

The first fourteen are charge- independent, 

0^=^'^^ = [1, a,, ■ a„ S,,, L ■ S, L^, LV^ • a,-, (L ■ S)^] ® [1, r, ■ r,] , (2.5) 



and the last four. 



Of-^^'^^ = [1, a, ■ a„ S,.j\ ® r,, + r,,- , (2.6) 



break charge independence. We will refer to the potential from the p = 15 — 17 terms as 
v*-^^ and from the p = 18 term as v^^^ . We note that in the context of isospin symmetry 
the CI, CSB and CD terms are respectively isoscalar, isovector and isotensor. 

The two-nucleon potential is supplemented by a three-nucleon interaction from the Ur- 
bana series of models 0, including both long-range two-pion exchange and a short-range 
phenomenological component: 

V^,k = + V,% . (2.7) 

The two-pion-exchange term can be expressed simply as a sum of anticommutator and 
commutator terms, 

y'^k = E + ySf ■ (2-8) 

cyclic 

Here 

V^^ = A,AX:„ XJ,}{n ■ r,, r, ■ n} , (2.9) 

and 
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yt^M = 4^2.[Xr„Xj;][r, ■ r,,r, ■ n] , (2.10) 

with X^j = Y{rij)(Ti -aj + T{rij)Sij as the basic one-pion exchange operator. The V^fj^ has no 
spin-isospin dependence, and in the Urbana model IX [^] its strength and A2tt are adjusted 
to reproduce the binding energy of and give a reasonable saturation density in nuclear 
matter when used with Argonne fig. 

The CD and CSB terms in H are fairly weak, so we can treat them conveniently as a 
first-order perturbation and use a wave function of good isospin, which is significantly more 
compact. Also, direct GFMC calculations with the spin-dependent terms that involve the 
square of the momentum operator can have large statistical fluctuations, as discussed in 
Ref. [|l|. Thus we construct the GFMC propagator with a simpler isoscalar Hamiltonian, 

^' = E^f' + E<+ E yi^k. (2.11) 

i i<j i<j<k 

where f is defined as 

'J 

<•= E 4{n^)0f, + v'cir,,) . (2.12) 

p=l,8 

The interaction v'^j has only eight operator terms, with operators [1, Ui ■ aj, S'jj, L ■ S] ® 
[l,Ti ■ Tj], chosen such that it equals the isoscalar part of the full interaction in all S and 
P waves as well as in the ^Di wave and its coupling to the ^5*1. The isoscalar part of 
the pp Coulomb interaction, t;^, is also included in H'. Detailed expressions are given in 
Ref. The v'^j is a little more attractive than Vij, so we compensate by using a V/jg. 
that is adjusted to keep {H') ^ {H); this should help prevent the GFMC propagation from 
producing excessively large densities due to overbinding. The small contribution of {H — H') 
is calculated perturbatively. 



III. VARIATIONAL MONTE CARLO 

The variational method can be used to obtain approximate solutions to the many-body 
Schrodinger equation, if\E' = i?\E', for a wide range of nuclear systems, including few-body 
nuclei, light closed-shell nuclei, nuclear matter, and neutron stars 0. A suitably parame- 
terized wave function, is used to calculate an upper bound to the exact ground-state 
energy. 



Ey = 





\H\^v) 







The parameters in \&v' are varied to minimize Ey, and the lowest value is taken as the 
approximate ground-state energy. 

Upper bounds to excited states are also obtainable, either from standard VMC calcu- 
lations if they have different quantum numbers from the ground state, or from small-basis 
diagonalizations if they have the same quantum numbers. The corresponding \l/y can then 
be used to calculate other properties, such as electromagnetic form factors and spectro- 
scopic factors P], or it can be used as the starting point for a Green's function Monte Carlo 
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calculation. In this section we first describe our ansatz for for the A = 8 nuclei and then 
briefiy review how the expectation value, Eq. (|3.1|) , is evaluated and the parameters of 
are fixed. 



A. Wave Function 



Our best variational wave function for the nuclei studied here has the form 



(3.2) 



i<j<k 



where the pair wave function, \E'p, is given by 



i<j 



(3.3) 



The Uij, U^j^ , Uijk, and U™^ are noncommuting two- and three- nucleon correlation oper- 
ators, and the 5 is a symmetrization operator. The Uij includes spin, isospin, and tensor 
operators, while the U^j^ has spin-orbit operators, refiecting the operator structure of the 
two- nucleon interaction, Eq. ( p.4|) . The Uijk is a nontrivial operator in spin-isospin space also 
induced by Vij, while the U™^ refiects the structure of the three-nucleon interaction. All 



these correlations are discussed fully in Refs. The two-body correlations are gener- 

ated by the solution of coupled differential equations with embedded variational parameters. 
We have found that the parameters optimized for the a-particle are near optimal for use in 
the light p-shell nuclei. Likewise, the best parameters for the three-body correlations are 
remarkably constant for different s- and p-shell nuclei, so they have not been changed from 
the previous work. 

The form of the totally antisymmetric Jastrow wave function, j, depends on the nuclear 
state under investigation. For s-shell nuclei we use the simple form 



i<j<k i<j 



(3.4) 



Here fciTij) and f^^j. are central two- and three-body correlation functions and for the a- 
particle, 



|<l>4(0000)) =A\p]p[n]n[) 



(3.5) 



The Jastrow wave function for the light p-shell nuclei is significantly more complicated due 
to the requirements of antisymmetry. Expressions for A = 6, 7 nuclei are given in Ref. [|l|; 
the present j is a straightforward extension. We use LS coupling to obtain the desired 
JM value of a given state, as suggested by standard shell- model studies |Tl|]. We also 
need to specify the spatial symmetry [n] of the angular momentum coupling of four p-shell 
nucleons [|T^. Different possible LS'[n] combinations lead to multiple components in the 
Jastrow wave function. We allow for the possibility that the central correlations fc{fij) 
could depend upon the shells (s or p) occupied by the particles and on the LS[n] coupling. 
The Jastrow wave function is taken as 
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I^j)='4{ n fijk n fssirij)Y[fsp{rk5)fsp{rm)fsp{rk7)f spirts) 

i<j<k<A i<j<4. k<4: 

1234:5678 )}• (3.6) 

LS[n] 5</<m<8 

The operator A indicates an antisymmetric sum over all possible partitions of the eight 
particles into 4 s-shell and 4 p-shell ones. For the two-body correlations we use /ss(r) = fc{r) 
from the ^He wave function, while 

fsp{r) = fc{r) + Csp{l - exp[-(r/4p)']) , (3.7) 

(3.8) 

and 

fjgir) = Ur) + cW(l - eM-{r/d,,f]) , (3.9) 

where we have supplemented the fc{r) with a Ion g-range tail. The Cgp-, dsp-} etc., are varia- 
tional parameters, whose values are given in Table |. For the three-body correlations, our 
best present trial function has the /^^-^ acting only within the s-shell. 

The LS[n] components of the single-particle wave function are given by: 

|$8(L5MJMTT3)i234:5678) = | $4(OOOO)i2340f (i?a5)0p''(i?a6)0p^(i?a7)0j^(i?a8) 

{ [^Imi (^a5)^lmj (^qg)^!™;' (^a7)^lmj" (^a8)]LMi 

X [X5 (|"^s ) X6 ( i"^'J X7 ( i ) X8 ( i m';') ] SMs } 

xh(it3)z/6(|4)^7(it:;)z.s(K)]TT3) • (3.10) 

The (f)p''^{Rak) are p-wave solutions of a particle of reduced mass ^rriN in an effective a-N 
potential: 

VaN{r) = V^^'{r) + V^^ir) . (3.11) 

They are functions of the distance between the center of mass of the a core (which contains 
particles 1-4 in this partition) and nucleon k, and again may be different for different LS 
components. For each state considered in the present work, we have used bound-state 
asymptotic conditions for the (f)p^ , even if the state is particle unstable. The Woods-Saxon 
potential 

V^^Hr) = Vp^'il + expC-^)r , (3.12) 

has variational parameters Vj^^, Rp, and ap whose values are given in Table |I[ The Coulomb 
potential is obtained by folding over nuclear form factors: 

r 2 1 — y 

-^expi-Xp)[2 + Xp + ^-^][l-y']-'} . (3.13) 
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Here Xa = \/Vlr/ra, Xp = \/Vlr/rp, and y = Va/vp, with the charge radii = 1.65 fm and 
Tp = 0.81 fm. This additional potential term can be used with strength Q = 0, |, or | for 
^He, ^Li, or ^Be, respectively, corresponding to the average Coulomb interaction between 
the a core and a p-shell nucleon. The wave function is translationally invariant, hence there 
is no spurious center of mass motion. 



The experimental spectra for A = 8 nuclei are shown in Fig. ||. The ground state 
of ^He is strong stable, but decays by (3~ emission with a half life of 119 ms. One excited 
state is identified at 3.59 ± 0.05 MeV [0], above the threshold for decay to ^He+2n. In the 



shell model, the {J'^;T) = (0+;2) ground state is predominantly a '^^~^^L[n] =^S[22] state, 
where we use spectroscopic notation to denote the total L and 5* and the Young pattern [n] 
to indicate the spatial symmetry. The (2+;2) first excited state is predominantly a ^D[22] 
state. We also allow for possible ^P[211] admixtures, which are the only other available 
configurations in the p-shell, in our After other parameters in the trial function have 
been optimized, we make a series of calculations in which the f^LS[n] may be different in 
the left- and right-hand-side wave functions to obtain the diagonal and off-diagonal matrix 
elements of the Hamiltonian and the corresponding normalizations and overlaps. We diago- 
nalize the resulting matrices to find the PlsIu] eigenvectors. The shell-model wave functions 
are orthonormal, but the correlated cire not. Hence the diagonalizations use generalized 
eigenvalue routines including overlap matrices. We also calculate the position of the three 
predominantly 3p[211] states, with ( J^; T) = (2+; 2), (1+; 2), and (0+; 2); none of these have 
been identified experimentally. The normalized Pisin] for these different states are given in 
Table |n|. 

The ground state of *Li is a strong stable (2"'"; 1) state that decays by (3~ emission with 
a half life of 838 ms, and is predominantly ^P[31] in character. The (l"*"; 1) first excited 
state at 0.98 MeV excitation is also strong stable with a 12 fs 7-decay, and is primarily a 
mix of ^'^P[31] configurations. The (3"'"; 1) second excited state at 2.26 MeV is just above 
the threshold for breakup into ^Li+n and is fairly narrow with a width of 33 keV, decaying 
by both 7 and n emission. A number of higher states have been identified, most having 
fairly large widths. Two special cases are the (4"*"; 1) stretched state at 6.53 MeV excitation 
with a width of 35 keV, which can only come from the ^F[31] p-shell configuration, and 
the (0^; 2) isobaric analog of ^He that occurs at 10.82 MeV with a width less than 12 keV. 
The possible p-shell components in the T = 1 states include the ^'^P[31], ^'^D[31], ^''^F[31], 
^S[22], ^D[22], and ^'^'^P[211] configurations. We include all but the lowest-symmetry [211] 
components in our and calculate all possible first, second and some third states of given 
( J'^; T) through the diagonahzation procedure discussed above. Table |T| gives a summary 
of the Pisin] amplitudes. 

The (O"*"; 0) ground state of ®Be is 92 keV above the threshold for breakup into two 
a-particles, with a width of only 7 eV; it is an almost pure ^S[4] configuration. The first 
and second excited states are very broad: a (2+;0) "^D[4] state at 3.04 MeV and a (4+;0) 
^G[4] state at 11.4 MeV — the spacing of an almost rigid rotor. Indeed, as discussed below, 
'^Be appears to have the intrinsic deformation of a 2a molecule. Higher in the spectrum 
is the famous pair of isospin-mixed (2+; + 1) states at 16.63 and 16.92 MeV which have 
widths of order 100 keV; the T = 1 component is the isobaric analog of the ®Li ground 
state. There are similar (l"*"; + 1) and (3~^; + 1) pairs near 18 and 19 MeV which have 
widths less than 300 keV. Many additional states have been identified above 18 MeV, most 
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of them with large widths, up to the (0"'"; 2) isobaric analog of ^He at 27.49 MeV. There are 
also some negative-parity states in this region, which we have not attempted to calculate. 
In constructing we use all p-shell states of symmetry [4] and [31], but neglect those of 



symmetry [22] and [211], as shown in Table |V 



The known experimental spectrum for is similar to ^Li, except that with the extra 
Coulomb repulsion, the (2+; 1) ground state is just barely strong stable, decaying by P+ 
emission with a half life of 770 ms. The (1"^; 1) first excited state is seen as a narrow 
resonance in ^Be+p scattering, and the (3^; 1) second excited state is much broader than 
its ®Li analog. The only other observed state is the (0"'"; 2) isobaric analog to ^He at 10.62 
MeV. The ground state for ®C is unstable against several possible breakup channels, having 
a width of ~230 keV, but its mass excess is known within 20 keV. We have calculated the 
energies of these states in ^B and so we can study the energy differences in the T = 1 
and 2 isomultiplets. 

The full A = 8 wave function is constructed by acting on the Eq. with the 



Uij, U^j^ , Uijk, and U™^ correlations. Because of the tensor and other correlations in Uij, 
many additional symmetry components, beyond the explicit p-shell states discussed above, 
are built up in the wave function. In principle, the Uij could be generalized to be different 
according to whether particles i and j are in the s- or p-shell, but this would require a larger 
sum over the different partitions and would increase the computational cost by an order of 
magnitude. 

For input to the GFMC algorithm, it is more efficient to use the somewhat simplified 
trial function 

i^t) = 5 n (i + 1/., + e ) 1^^) > (3-14) 

where U^^^ is a truncated three-nucleon interaction correlation based on the short-range 
V^i^ term and on the anticommutator part of the two-pion exchange, ^j^"^, which can be 
reduced to operators that depend only on the spins and isospins of nucleons i and j. Thus 
the sum over k can be made, leaving a two-body spin-isospin operator that can be combined 
with Uij] the result is calculable with only a little more effort than just Uj alone. This trial 
function gets the bulk of the energy, as shown below, but for about half the computational 
effort of the full 



B. Energy Evaluation 



The energy expectation value of Eq. ( |3.1|) is evaluated using Monte Carlo integration 
A detailed technical description of the methods used here can be found in Refs. ||10|,|15|,|16 



Monte Carlo sampling is done both in configuration space and in the discrete order of 
operators in the symmetrized product of the pair wave function by following a Metropolis 
random walk. The expectation value for an operator O with the full wave function \l/yis 
given by 



(3.15) 
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where we have introduced a probabihty distribution, Wpq(Tl), based on the approximate 
wave function \l/p of Eq. ( |3.3|) , 



W„{K) = |Re((v^l,,,(R)vl/p,,(R)))| . (3.16) 

The subscripts p and q specify the order of operators on the left and right hand side of the pair 
wave functions, while the integration runs over the particle coordinates R = (ri, r2, . . . , r^). 
This probability distribution is much less expensive to compute than one using the full wave 
function of Eq. { \i.'2\ ) with its spin-orbit and operator-dependent three-body correlations, but 
the denominator of Eq. (|3.15|) is typically within 1-2% of unity. Expectation values have 
a statistical error which can be estimated by the standard deviation a in either gaussian 
approximation or by using block averaging schemes. 

Our wave functions are vectors of 2^ x I{A, T) complex numbers, 

vI/(R) = ^^„(R)|a) , (3.17) 

a 

where the are the coefficients of each possible spin-isospin state |a) with specific third 

components of the spins of each nucleon and the desired total isospin. The CD and CSB force 
components are sufficiently small that we do not worry about isospin mixing in the wave 
function, and the expression for the number of isospin states, I{A,T), is given in Ref. [|l|. 
This gives arrays with 3584, 5120, and 7168 elements for ^Be, ®He, and ^Li, respectively. The 
spin, isospin, and tensor operators Of~'^'^ contained in the two-body correlation operator 
Uij, and in the Hamiltonian are sparse matrices in this basis. 

Expectation values of the kinetic energy and spin-orbit potential require the computation 
of first derivatives and diagonal second derivatives of the wave function. These are obtained 
by evaluating the wave function at 6A slightly shifted positions of the coordinates R and 
taking finite differences, as discussed in Ref. |T^. Potential terms quadratic in L require 



mixed second derivatives, which can be obtained by additional wave function evaluations 
and finite differences. A rotation trick can be used to reduce the number of additional 



locations at which the wave function must be evaluated [jT7|. 

In addition to calculating energies, we evaluate ( J^) and (Jz) expectation values to verify 
that our wave functions truly have the specified quantum numbers. Another check is made 
on the antisymmetry of the Jastrow wave function by evaluating, at an initial randomized 
position. 



(3.18) 



where P^j"'^ are the space, spin, and isospin exchange operators. This value should be 
exactly zero for an antisymmetric wave function, and it is in fact less than 10~^ for each 
pair of particles in each nuclear state that we study. 

A major problem arises in minimizing the variational energy for p-shell nuclei using the 
above wave functions: there is no variational minimum that gives reasonable rms radii. For 
example, the variational energy for ^Li is slightly more bound than for ^He, but is not more 
bound than for separated ^He and nuclei, so the wave function is not stable against 
breakup into a + d subclusters. Consequently, the energy can be lowered toward the sum 
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of ^He and energies by making the wave function more and more diffuse. Such a diffuse 
wave function would not be useful for computing other nuclear properties, or as a starting 
point for the GFMC calculation. 

On the basis of our work in A = 6-7 nuclei, we believe that part of the problem is a 
fault of the Hamiltonian, but most of it is due to the presence of small admixtures of highly 
excited states in the trial function 0]. In that work, we constrained our search for optimal 
variational parameters by requiring the resulting point proton rms radius, Vp, to be close to 
the experimental values for ^Li and ''Li ground states. Then we allowed only small variations 
in the construction of the ^He and ^Be ground states and all the excited or resonant states, 
for which there are no experimental measurements of the charge radii. 

In the present work, we have the additional complication that all the A = 8 nuclei are 
sufficiently short-lived that precise experimental determinations of their radii exist do not 
exist. Consequently, the variational parameters in Table | are chosen very close to those 
of our previous work, with only a systematic reduction in the depth of the Woods-Saxon 
potential well as L increases, and an increase in the tail of the f^{r) correlation as the 
spatial symmetry declines, as shown in Table |. These gradual changes help to insure that 
the radii of excited states increase as the excitation energy increases. 

The last step is always the diagonalization to determine the PlsIh] niixing coefficients 
of Tables 0-|IV|. Shell model lore tells us that the lowest state of any given ( J'^; T) will 
be the state with maximal spatial symmetry and smallest L that can be formed from the 
allowed couplings, e.g., the ^S[4] ground state in ^Be or the ^P[31] ground state in *Li. 
For the purposes of obtaining a variational upper bound and a GFMC starting point, we 
could settle for a constructed using only that LS[n] component. However, by using 
more components, we can gain a significant amount of energy in some cases and this gain 
persists in our GFMC propagations. For the ^He ground state, when the dominant ^S[22] 
piece is supplemented by the ^P[211] term, the energy gain is 1.2 MeV. In the case of ^Li 
ground state, where the dominant term is '^P[31], addition of the three other contributing 
[31] symmetry states lowers the energy 0.8 MeV; addition of the one [22] state (the next 
highest spatial symmetry) gives only an additional 0.1 MeV. However, for the ^Be ground 
state, with the dominant ^S[4] term, addition of the one [31] symmetry piece gives only 0.1 
MeV. For this reason, we feel reasonably confident in truncating the '^Li and '^Be p-shell 
bases after the top two symmetry states. 

IV. GREEN'S FUNCTION MONTE CARLO 

A detailed description of the nuclear GFMC method, and many tests of its accuracy, are 
given in Ref. |l|. In this section we present a brief review of the method and then describe 
two improvements that have been made since that publication. In most of this section, we 
will not make the distinction between H' and H; the reader should remember that in fact 
we use the simpler H' in our GFMC propagator and evaluate {H — H') perturbatively. 
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A. Review 



The GFMC method starts with the trial wave function, of Eq. ( |3.14| ), and projects 
out of it the exact lowest energy state with the same quantum numbers, \E'o- 



= lim "^{t) 



(4.1) 



T 



-{H-Eo)At 



(4.2) 



Here we have sliced the imaginary propagation time, r, into a number of small time steps. 
At = r/n. The small-time-step Green's function, Gq,^(R, R'), is a matrix function of R 
and R' in spin-isospin space, with matrix elements defined as 



G„^(R,R') = (R,a|e-(^-^«)^"|R',/5) . 
Then \E'(Rn, t) is given by 

^(R„, t) = J G(R„, R„_i) ■ ■ ■ ^(Ri, Ro)^T(Ro)c?Rn-i ■ ■ ■ c/Ri^Ro . 

The small-time-step propagator used in Ref. [H is 



(4.3) 



(4.4) 



Ga^{R, R') = e^o^"Go(R, R') exp 



x(a|/3(R)|7)(7l 



\5){5\h{^m , 



(4.5) 



where gij{Yij,v[j) is the exact two-body propagator, go^ij{rij,r'^j) is the free two-body prop- 
agator, (^0(1^, R-') is the free many-body propagator, a, j3, 7, 6 are spin-isospin state indices, 
and summation over 7, S is implied. There is also an implicit sampling of the order of pairs 
in the symmetrized product of Eq. ( [4.51 ). The construction of the exact two-body propaga- 
tor is described in Ref. [|T|. The influence of the three- nucleon potential on the many-body 
propagator is broken into two pieces: the scalar which is easily exponentiated, and the 
V^l which is a more complicated operator in spin-isospin space. The simplest treatment of 
this term in the TNI is to expand to first order in At, 



h(R) 



1 



At 



2-^R) 



(4.6) 



in fact we use a more efficient procedure described below in subsection E. 

The integrals in Eq. ( [4.4| ) are evaluated stochastically by averaging over a set of n-step 
paths, P„ = Ro, Ri, ■ ■ ■ , Rn. The paths are chosen by first sampling a set of positions, 
Ro, using a probabihty function based on \E't(Ro)5 and then sequentially sampling the free 
Green's functions, Go(R'j+i5 R-i), to generate Ri+i from Rj. We thus obtain a spin-isospin 
vector ^E'(Pn) which is one sample of the integrand of Eq. ( [4.4| ). A GFMC "configuration" 
consists of the position, R„, and the vector \E'(P„). Branching and importance sampling. 
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described in detail in Ref. are used to obtain samples with probability proportional to 
the scalar importance function /: 



/[^(P„),^T,p(Rn)] = |ReE^„(P„)VT,p,„(Rn)]|+eEl[V'"(P-)Vr,p,a(Rn)]| , (4.7) 

a a 

where the sums run over the spin-isospin states a, and \1/t,p is the trial wave function 
evaluated with a specific choice of pair operator orders p. Here e is a small constant (~ 0.01) 
that ensures a positive-definite importance function so that diffusion can take place across 
nodal surfaces. 

The GFMC method allows one to compute "mixed" expectation values: 



(0)M^.ed - ^^^^^ ' (4-8) 



/rfR„^URn)^(Rn,r) 



(4.9) 



Because H commutes with the propagator, {H{t)) Mixed is an upper bound to Eq and ap- 
proaches Eq from above. However expectation values of operators that do not commute 
with H are extrapolated using 

^ (^(r)|0|^(r)) 
(^(r)l^(r) ) ' 

~ {0{r))M^xed + \{0{r))^uxed ' {0)t\ , (4.10) 

where 

(0>,= <#21M. (4.11) 



In the following we address refinements which have been made to the GFMC algorithm 
since Ref. |jl| . These refinements are important in order to make calculations of larger nuclei 
feasible. 



B. Constrained Path Algorithm 

Diffusion or Green's function Monte Carlo simulations of many-fermion systems generally 
suffer from the so-called "fermion-sign problem" . In essence this results from stochastically 
evaluating matrix elements of the form encountered in Eq. ( [4.9| ). The Monte Carlo techniques 
used to calculate the path integrals leading to \l/(R„,r) involve only local properties, while 
antisymmetry is a global property. This leads to integrands in Eq. ( |4.9| ) that have oscillating 
signs at large r, which cause the statistical error to grow exponentially with imaginary 
time. The problem is not insurmountable for light systems, because the states are fairly 
well separated and one can propagate for a substantial imaginary time without having 
unacceptable statistical errors. However, the sign problem also grows exponentially with 
particle number, as the interchange of any pair of nucleons causes a change in sign for the 
matrix element. 
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The sign problem significantly limits the maximum r that we can use in the simulations; 
hence we invoke an approximate technique to deal with it. The approximation involves 
keeping only a subset of the paths in evaluating the integrals, and using the knowledge 
gained in the VMC calculations to choose the subset of paths. It is closely related to 
methods used previously in condensed matter and elsewhere; they generally go by the name 



of constrained-path techniques 0,|T^. Some details of the algorithm, however, are special 
to the nuclear physics case. 

The basic idea of the constrained-path method is to discard those configurations that, 
in future generations, will contribute only noise to expectation values. If we knew the exact 
ground state l^^o), we could discard any configuration for which: 

^(P„)^^o(Rn)=0, (4.12) 

where a sum over spin-isospin states is implied. The sum of these discarded configurations 
can be written as a state I^E'd), which obviously has zero overlap with the ground state 

(^,|^o) = . (4.13) 

The contains only excited states and should decay away as r — oo, thus discarding it 
is justified. However, in general the ground state \E'o is not known exactly and hence the 
constraint is imposed approximately using \E^t in the place of \E'o. 

In Green's function or Auxiliary Field Monte Carlo (AFMC) p0|-p2|, the overlap of the 



configuration with the trial state evolves smoothly with time. The change of the configu- 
rations ^(Pri) per time step scales with y/ At, which can be made arbitrarily small. If the 
wave function of the system is a purely real scalar quantity, any configuration which yields 
a negative overlap must first pass through a point at which and hence the overlap is 
zero. Discarding configurations at this point is sufficient to stabilize the simulation and pro- 
duce an approximate solution \I'c to the many-fermion problem. It solves the many-body 
Schrodinger equation with the boundary conditions imposed by the nodes of "^t, and is 
known as the fixed-node approximation [p^-p6| . The approximate \E'c is the best wave func- 
tion (in the sense of lowest energy) with the same nodes as "^t- The discarded configurations 
are orthogonal not only to the trial state but also to the solution of the constrained problem 
\E'c and it has been shown that this method produces an upper bound to the ground state 
energy 

More generally, and particularly in nuclei, the trial wave function is a vector in 
spin-isospin space, and there are no coordinates for which all the spin-isospin amplitudes 
are zero. This is also true in AFMC, where the propagated configurations describe a fully 
antisymmetric product of single-particle wave functions. In both these cases, the evolution 
of the wave function remains smooth. In AFMC it is still possible to use a constraint of the 
form given in Eq. ( [4.13| ). However, the discarded configurations, while orthogonal to the 
trial state \E't, are not necessarily orthogonal to the the solution in the constrained space 
"ifc- For this reason, the method does not produce an upper bound to the true ground-state 
energy [0. 

A new difficulty specific to the nuclear problem is that the overlap \E'T,p(R'n)^^(Pn) is 
complex, and is estimated stochastically with a randomly selected order, denoted by the 
subscript p, of the operators in "^t, Eq- ( |3.14|) . These overlaps do not evolve smoothly and 
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pass through zero. Therefore we can satisfy the constraint Eq. ([4.13|) only on the sum of 
discarded configurations {"^d), but not for individual configurations as in Eq. ( 4.12|) . 



The fluctuations in '^T,p, due to the sampling of the pair operator orders, are small 
compared to the wave function itself, so we define an algorithm for discarding configurations 
which resembles as much as possible the fixed-node or constrained-path algorithms described 
above. Configurations in the GFMC are obtained with probability proportional to the 
importance function It,p = /[\1/(P„), \1/t,p], Eq. ([4.7]), which depends upon p. For each of 
them the overlap Ot,p is defined as: 

Ot,p = Re[*(PO^^T,p(Rn)] , (4.14) 

so that the required constraint is the sum over discarded configurations of Ot,p/It,p = 0. 
We define a probability P[\I^(P„), \l/r,p(Rn)] for discarding a configuration in terms of the 
ratio Ot,p/It,p- 

P[^(P„),^T,p(Rn)] = O/I > a, 

= ^ «e > O/I > (4.15) 

= 1 O/I < p,. 

According to this algorithm configurations with O/I less than f3c are always discarded, con- 
figurations with O/I greater than Oc are never discarded, and there is a linear interpolation 
in between. The values of etc and Pc are held constant during the simulation. At each step, a 
configuration is discarded with a probability that depends upon the order of pair operators 
p. Hence an unbiased estimate of the overlap is required. It is obtained by choosing 

independent samples p' of the pair operator order to compute the overlap: 

i^dl^T) = E Ot,p'/It,p' . (4.16) 

discarded 
configurations 

One or two preliminary runs are sufficient to adjust the constants and Pc such that this 
overlap is zero within statistical errors. When this condition is satisfied, it is trivial to show 
that the growth estimate of the energy, obtained from the growth or decay of the population 
with imaginary time, is identical to the mixed estimate. These two energy estimators are 
described in detail in Ref. 

Note that in absence of the term with e in the importance function I^p, Eq. ([4.7|) , the 
ratio Ot,p/It,p is ±1. In constrained path calculations, a significantly larger value of e, 
typically 0.15, is used to reduce the fluctuations in Ot,p/It,p- For light nuclei, we have tried 
a variety of Uc and Pc adjusted to the constraint without finding any statistically significant 
differences in the results. Indeed, even preliminary estimates of ac and Pc which do not 
precisely satisfy the constraint condition Eq. ( |4.13| ) yield results indistinguishable from the 
results with finely tuned constraint parameters. Typical values of ac and Pc for A = % are 
0.1 and 0.2 respectively. 

In principle, two parameters are not needed to adjust the constraint to yield zero overlap. 
We were motivated to try this approach by the fact that, in standard fixed-node calculations, 
one defines a propagator which goes exactly to zero at the node. This results in some con- 
figurations being discarded which have small positive overlaps, as well as all configurations 
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which have negative overlaps. In the fixed node scheme, this produces a result of higher 
order in the time step than simply discarding configurations with negative overlaps. For this 
reason both and /?c are chosen to be greater than zero. The algorithm described above is 



related to one proposed by Sorella in the context of condensed matter simulations l^g], but 
simpler in that it uses only the overlap and the constraint parameters are held fixed during 
the simulation. 

As discussed in the next subsection this method yields results with stable statistical 
errors independent of r. For A < 7, the calculated energy is very close to that obtained 
without constraints in most cases; however in some cases, particularly those with poor ^l/r's, 
the calculated energy can be significantly below or above the true unconstrained result. 
The amplitudes of the excited states with energies E'j, in decay exponentially as 
exp(— — Eo)t) in the unconstrained ^I/(t) without any change in their phase. Therefore 
these give a positive contribution to the mixed energy, making it >Eq. The equality is 
obtained in the limit r— >-0. The constrained energy can have additional errors if, due to a 
poor choice of high-energy excitations get reintroduced in \l'(r) — Such excitations do 
not necessarily have the same phase as in and can give contributions to the constrained 
E{t) of either sign. Such errors can be easily detected by propagating \l/(r) — without 
constraint for a limited number of steps n^. Most of the dependence of the calculated 
energy is eliminated by fairly small n„ without a significant increase in the statistical error. 
Thus the wave function \l/(7(r, n^) used to estimate energies and other observables in the 
present calculations, is obtained from: 

= nAr,ri„) = {exp[-ffn„Ar]}„{exp[-iJ(n - n„)Ar]}e^T , (4.17) 

where {■ ■ ■}c signifies propagation with the constraint and {• ■ •}„ indicates normal propaga- 
tion without constraints. As is shown in the next subsection, the simulations are stable for 
arbitrary n, and a fairly small n„ {huAt < 0.01 MeV~^) is sufficient to eliminate the \I/t 
dependence. Evaluation of matrix elements of this wave function can be implemented very 
easily on the computer by, after each propagation step, labeling those configurations which 
are to be discarded and then retaining them in the simulation for Uu more steps. In principle 
this could be used to evaluate the overlap of the discarded configurations with after the 
unconstrained propagation. This could result in better values of the constraint parameters, 
however we find no change in the overlap for the values of we have considered. 



C. Tests of Constrained Path 

We have tested the constrained path algorithm in a variety of light nuclear systems, 
studying the dependence upon constraining wave function and also the convergence of the 
results obtained by relaxing the constraint. In some cases bad trial wave functions were used 
to test the algorithm under extreme conditions. In this subsection we present results for the 
^Li and ^He ground states and for 8 neutrons bound in an external well. The tests for ^Li 
and 8 neutrons were made using H = H', Eq. ( p.ll|) , with no three-body potential. This 



eliminates uncertainties from the extrapolation oi H — H', and allows us to have just the 
exact two-body propagator, gij. Eliminating the three-body potential also results in faster 
calculations allowing smaller statistical errors to be achieved. In the following all energies 
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are in MeV and imaginary times in MeV~^; for simplicity the units are omitted in most 
cases. 

Figure ^ shows the ground state energy of ^Li calculated with two choices of 
[Eq. ( p.l4|) ]: (1) which is the full wave function with no Ujj^^ because there is no 
Vijki and (2) \I^a-, a simplified wave function obtained by removing the tensor pair correla- 
tions leaving only spin-spin and isospin correlations. The subscript t on emphasizes that 
it has the essential tensor correlations. In this and the following two figures, constrained- 
propagation results with zero and finite ?t,„ are shown as open symbols while solid symbols 
show completely unconstrained results or unconstrained continuations of constrained prop- 
agations. Averages of the last few energies are shown by solid lines with the corresponding 
statistical errors indicated by the surrounding dashed lines; the length of the line indicates 
the range of r included in the average. Consider first the results using the full (\l/f) wave 
function shown in the upper part of Fig. § these are typical of our production calculations. 
The solid circles show an unconstrained propagation from r = to r = 0.1 of 200,000 config- 
urations; the statistical errors on the last points are growing rapidly, but the energy seems to 
be fairly independent of r beyond 0.05. The average from r = 0.06 to r = 0.1 is -28.16(12), 
as shown by horizontal solid and dashed lines. It should be close to the correct answer for 
this Hamiltonian. The open triangles show constrained results with n„ = for the same 
case. This simulation is stable with much smaller statistical errors than the unconstrained 
calculation and could be extended to an arbitrarily large r; here we find no significant change 
after r = 0.05. The average over 0.14 < r < 0.20 is -27.94(2). We released the constraint at 
r = 0.20 and continued the propagation to r = 0.27, as shown by solid triangles. The errors 
grow, and the average over this entire region is -28.15(5). It is likely that this result, which 
is consistent with the unconstrained result of -28.16(12) but with smaller statistical error, is 
the most accurate available. A constrained propagation with = 20 is shown by the open 
squares; only 50,000 configurations were used. This simulation is very similar to the = 
one; the average is -28.01(5). These results suggest that constrained propagation with our 
typical \l/f, using n„ = 20, may lead to an error of +0.14(7) (or ~0.5%) in the binding energy 
of '^Li. 

The results in the lower part of Fig. |^ were obtained using the ^ ^ trial wave function. 
This is an extremely bad ^j- because, with no tensor correlations, this wave function results 
in an identically zero expectation value for the tensor potentials. Thus while the r = 0.0 
expectation value for '^t is -23.95(4), that for is +31.1(1). The solid circles show an un- 
constrained propagation of 400,000 configurations for this case; the GFMC has managed 
to improve the r = 0.0 energy by 59 MeV to -28.0(2), but the statistical errors are large 
even after using twice the number of configurations as for the unconstrained propagation 
with The open triangles show = constrained propagation of 200,000 configurations. 
This propagation becomes stable after r = 0.12, but it has converged to the significantly 
overbound result of -30.2(1). The solid triangles show the results of releasing the constraint 
at r = 0.20; the energy immediately goes up to the correct value with an average of -28.0(2). 
Finally the open squares show n„ = 20 constrained propagation of 400,000 configurations. 
The repeated use of 20 unconstrained steps gives a stable propagation that is quite accurate; 
the average is -28.27(7) with an error of -0.12(9). The unconstrained continuation of this 
calculation, shown by the solid squares, makes no significant change. These results show the 
need to have unconstrained steps before evaluation of the energy or any other observable. 
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The second test case is a system of 8 neutrons bound in an external one-body potential. 
We have previously reported results for such systems as a basis for comparing Skyrme models 



with microscopic calculations based on realistic interactions . As in the previous example, 
the neutrons interact via the v'^j, with no three- nucleon interaction. Because systems of 
neutrons are not self binding, the Hamiltonian also includes an external one-body potential 
of Woods-Saxon form, 

the parameters are Vq = —20 MeV, tq = 3.0 fm, and = 0.65 fm. Neither the external well 
nor the internal f potential are individually attractive enough to produce a bound state of 
eight neutrons, however the combination does produce binding. 
The Jastrow for neutron drops is just 

\^J) = l[fcir^J)\^AiJM)) . (4.19) 

i<j 

where $a is a Slater determinant of single-particle orbitals. We compare results obtained 
with the full ^'t, again refered to as because of the tensor operators, shown in the upper 
part of Fig. ^ with those obtained with the shown in the lower part. 

The solid circles show results of unconstrained propagation using 400,000 configurations 
of the and 560,000 configurations of the trial wave functions, respectively. The errors 
are growing rapidly by r = 0.06, and the calculations do not appear to have converged. 
Averages of the last four energies are -39.3(1) and -38.9(1) for \E'f and \E'j respectively, 
but these are clearly just upper bounds. Constrained propagations using the trial wave 
function and = and = 20 are shown by the open triangles and squares, respectively 
in the lower part of the figure. These have stabilized beyond r = 0.1, with end averages of 
respectively -39.58(4) and -39.79(6). Both propagations were continued without constraint; 
their results are shown by solid triangles and squares. The energy remains stable but the 
errors grow rapidly; the averages are -39.86(11) and -39.85(34). In this case constrained 
propagation with n„ = 20 gives results that are below the best unconstrained (starting 
from r = 0) upper bound with by 0.5(1). The constrained result is consistent with 
unconstrained continuations beyond r = 0.25, indicating its reliability. 

A different situation obtains with the variationally better trial wave function; the 
r = 0, VMC energy for ^t, -35.30(4), is lower than the -30.62(9) given by ^j, but the 
difference is not as large as in the case of ^Li. Constrained propagation using and = 
is shown by the open triangles in the upper part; it stabilizes beyond r = 0.1, but at a 
value that is 1 MeV too high. The results of releasing the constraint are shown by the solid 
triangles; the energy immediately drops 1 MeV and the final average is -39.71(10). Finally 
constrained propagation for and = 20 is shown by the open squares. The results 
are stable with an average of -39.74(5), in agreement with the best results obtained in the 
previous paragraph. Releasing the constraint does not result in any significant change, as 
is shown by the solid squares. This case also confirms the need to release the constraint 
before measurements, with the surprising result that a variationally better wave function 
may not necessarily provide better constraints to guide the GFMC. Nevertheless, constrained 
propagation with = 20 gives results with presumably less than 1% error, for both trial 
wave functions. 
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Finally, in Fig. ^ we show results for the ^He ground state with the full Argonne vis plus 
Urbana IX three-nucleon interaction. We are presenting this case because we find that '^He is 
the most difficult nucleus we have studied, in terms of obtaining reliable error estimates, and 
the need for n„ > 0. All of the calculations are made with the full of Eq. (|3.14|) . In this 
figure the results of a standard calculation without constraint (solid circles) are compared 
to those with constrained-path = (open triangles) and n„ = 20 (open squares), and 
their unconstrained continuation beyond r = 0.2 shown by solid triangles and squares. 
The average unconstrained energy in the 0.03 < r < 0.06 range, -26.1(3), is clearly an 
unconverged upper bound. The n„ = and 20 end averages, -26.89(9) and -27.16(15), are 
below it by ~1 MeV. The averages of the results obtained after the constraint is released at 
r = 0.2, -27.5(4) and -26.9(2), are not significantly different from the constrained averages. 

From these, and many other, tests we conclude that constrained propagation, including 
unconstrained steps prior to measurement, yields results that are reliable. In the A = 8 
nuclei and neutron drops they are up to 4% below the unconverged upper bounds that 
can be obtained by unconstrained propagation up to r = 0.06, and have smaller statistical 
errors. The constrained path results with riu ~ 10 to 20 are stable within 1% with respect 
to reasonable, and in case of ^Li even unreasonable, changes in used to constrain the 
paths. Our present practice is to use the full Eq. ( |3.14[ ) in constrained path calculations 
with n„ = 20 for the neutron-rich He isotopes and neutron drops and n„ = 10 for all other 
nuclei where even = calculations seem to be fairly accurate. 

D. Resonance states in constrained-path algorithms 

Since the constrained-path algorithm is stable to large imaginary time, it is useful to 
consider the large-r behavior of the energy of resonant states. The bound-state simulations 
described above yield asymptotically stable energies out to very large imaginary times. 
Resonant states, however, are more delicate. In principle these states should decay to 
separated clusters, and in fact this does occur with the constrained-path algorithm. The 
rate of this decay presumably depends not only on the resonance energy but also on the 
width of the state, and the trial state wave function used to impose the constraint. 

We have studied several cases including the two unbound p-wave states in ^He. The 
GFMC energy consistently decreases with r for both the J'^ = l/2~ and the 3/2~ states. 
The 1/2^ state, lying higher in energy and with a larger width, decays more quickly than 
the 3/2~ state. We have verified that the system breaks apart into a separated a plus n by 
studying the pp pair distribution function, which depends only upon the internal structure 
of the a-particle. This distribution remains constant out to the largest imaginary times (3 
MeV~^) studied, while the pn and nn distributions steadily become broader. 

We have also studied three low-lying states in ^Li out to large r. The = 1+ ground 
state is stable, the 3^ excitation is a narrow resonance, while the 2"*" state lies higher in 
energy and is broader. The GFMC energies of these states are displayed in Fig. ^ as a 
function of r. The ground state is clearly stable out to the largest imaginary times. The 
narrow 3^ shows a plateau in E{t), decreasing only modestly in energy, though of course 
it will eventually also decay to the a+d threshold energy. The 2+ state is much broader 
and decays much more quickly in imaginary time. Therefore, comparison of the energies 
obtained for such broad states with experiment could be misleading. The a-|-d threshold 
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energy for this Hamiltonian is shown in the figure; it is clear that both excited states are far 
from convergence to this energy. 

It is possible to compute energies that may be directly compared with experiment for res- 



onant states that have only a single two-body channel available for breakup [Q. Scattering 
observables, including scattering length, effective range, and phase shifts, can be calcu- 
lated directly by imposing scattering boundary conditions on the asymptotic wave function. 
This is quite important for studying a variety of interesting low energy reactions, including 
parity-violation and important astrophysical reactions, but has not been considered in this 
work. 



E. Three-body Propagator 



The small-time-step propagator of Eq. (|4.5|) involves two complete evaluations of the 
three-body potential (two sums over all triples) and one product of all pair propagators for 
each time step. Thus as the number of nucleons is increased, the time spent in the three- 
body part of the propagation becomes a larger and larger fraction of the total time for the 
calculation. For this reason it is desirable to find a less costly treatment of the three-body 
propagator. 

It was noted in Eq. ( p. 141 ) that the full can be replaced with a Uf^^^ that omits 



the commutator part of V"J^ with very little degradation of the variational energy, and that 
the resulting correlation is an operator in the spins and isospins of only two nucleons. This 
led us to consider the combined two- and three-body propagator (we omit the spin-isospin 
indices) : 



G(R,R') = e^°^^5n 



i<j 



90; 



«JV^«J5 '■ij 



r') 

'■11 J 



(4.20) 



1 - ^ E «C/(R) 



9ij[^iji ^ 
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X exp 



(4.21) 



The V^^!^^ is defined in Eq. (|2.9| ) as the anticommutator part of VM, and a is chosen so that 



(4.22) 



typically a = 1.6. This gij can be reduced to a single 4x4 operator in i,j spin space 
for each of the two isospin states, and thus takes little more time to evaluate than just 
gij. Propagation with just G gives much the same results as with the much more costly 
propagator of Eq. ([4.5|), but we attempt to make a more reliable propagation by using the 
following sequence as the basic propagation step: 



G(Ri+n, R'i+n-l)G'(R'i+n-l; ) 



G(R.i+n/2+l, Rj+ri/2) 



x{l-nAr Y: [V^.5(R.+„/2)-«V;5f^^(R.+n/2)]} 



i<j<k 

xG{R. 



■i+n/2, Ri+n/2-l) 



' (-^(Rj+i, Hi) , 



(4.23) 
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where n is a small (typically 4) number of steps. Here we go from Rj to Rj+n by making 
n/2 steps with G; then applying n times the correction due to the difference of the complete 
V'^'^ and the approximate aV'^'^'^^ both computed at the position Rj+„/2; and finally making 
another n/2 steps to Rj+n- 

Table |V| shows the reliability of this method for ^Li(gs). The first line gives results with 
the old method, Eq. (|4.5| ), and the next three lines use the new method with increasing 
values oi n [n = oo means that no correction for the approximation in G is made with 
Eq. ( [4.23| )]. The last line shows results using just a two-body propagator; that is, the three- 
body potential is included only perturbatively. Clearly a three-body propagator is essential 
for obtaining results with < 1% error, however the G is reliable to better than 1%. 

F. Numerical Evaluations 

The p-shell GFMC calculations reported here were made with constrained-path prop- 
agation to r = 0.2 MeV~^ using 400 steps of Ar = 0.0005 MeV~^. Energies and other 
observables were evaluated every 20 steps and the last 7 (r > 0.14 MeV-^) values were 
averaged. At least 10 unconstrained steps were taken before the observables were computed 
(20 steps were used for *He). The A = 8 results used 10,000 to 20,000 initial configurations. 
Once the propagation has stabilized (typically by r = 0.1 MeV~^), the constraint removes 
from 0.7% [for *Be(gs)] to 1.7% [for ®Be(3"'")] of the configurations at each propagation step 
(the percentages for ^He and ^Li are in this range). The removed configurations are re- 
placed by new ones generated by branching with the average branching probability chosen 
to maintain an approximately constant population. 

Most of the GFMC calculations reported here were made on the 128-node SGI Origin 
2000 in the Mathematics and Computer Science division of Argonne National Laboratory. 
The individual nodes in this machine are 250 MHz RIOOOO processors which (when the 
processors are not being shared with other users and memory remains local to each processor) 
deliver sustained speeds of > 200 MFLOPS for the 8-nucleon calculations. This results in 
a 10,000 configuration calculation for ^Be taking about 460 node-hours. The approximate 
times required for other nuclei can be determined from Table |VI[ The columns of this table 
show the number of nucleons, the number of pairs, and the size of the spin-isospin vector. We 
find that the total calculational time is proportional to the product of these three numbers; 
this product (scaled to 1 for '^Be) is given in the last column. 

Typically we have computed the M = J state of a given nucleus; this allows us to directly 
evaluate the spectroscopic quadrupole moment and magnetic moment. Recently we have 
realized that if we compute the M = state of even-A nuclei, the size of the spin-isospin 
vector can be reduced by a factor of two. This is done by observing that for even A, 

vl/_,,,_,,,.,_,,(J,M = 0) = (-l)i-^+-^-^^^vl>;^^^ ,„^^JJ,M = 0) , (4.24) 

where Sj is the spin-projection of nucleon i and Ms = J2 Si- Thus only half of the spin- 
isospin vector needs to be computed and stored, resulting in a saving of half the computer 
time for even A > 8 nuclei. This saving is not included in the times discussed in the previous 
paragraph. Unfortunately, the corresponding relation in odd-A nuclei relates \E'_s^(J, M) to 
"^SiiJ, —M) and thus does not reduce the computational effort. 
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V. ENERGY RESULTS 



A. Ground-State and Excitation Energies 

In this Section we present VMC and GFMC results for the Hamiltonian consist- 
ing of Argonne Vig plus Urbana IX. Figure ^ shows the E{t) for the lowest T = 0, 
J'^ = 0"*", 2+, 4+, 1"^, 3"^ states of *Be. The solid and dashed lines show the average en- 
ergy and its statistical error; these numbers are reported in Tables |V11| and |V114 The E{t) 
are 

E{t) = {H'{r))Mi,ed + {H{r) - H'{r)) , (5.1) 

where the {H — H') is perturbatively extrapolated by Eq. ( |4.iq ). Figures and || show 
the corresponding '^Li and ®He results. In all cases for which the states are experimentally 
narrow or stable, the E{t) rapidly decreases for small r and stabilizes before r = 0.1 MeV~^. 
It is clear that E{t) has not converged by r = 0.2 MeV~^ for the experimentally very broad 
^Be(4"'") state, making it not possible to determine an accurate excitation energy for this 
state. The other broad states pBe(2~*') and ®He(2"'")], and the experimentally unknown 
^He(l"'") state seem reasonably converged. 

Table |VI1| shows the computed and experimental ground-state energies. The errors shown 
in parentheses are only the Monte Carlo statistical errors; the systematic errors discussed 
in the previous section and in Ref. |^ could add an additional 1% to the GFMC error. 
This paper is our first formal publication of A = 8 results. Results for A = Q and 7 were 
published in Ref. [D] ; most of the corresponding values in Table |V11| have been recomputed 



using improved VMC wave functions and GFMC propagation. Most of the VMC values have 
not changed significantly from those in Ref. |I| ; the exception is ^He which is now constrained 
to what we think is a more reasonable rms radius and actually has a higher energy. The 
A = Q and 7 GFMC energies have now all been computed with constrained propagation 
to r = 0.2 MeV^^ instead of the unconstrained propagation to only r = 0.06 MeV^^ used 
in Ref. [|l|. This has resulted in a lowering of ^'^He energies by 0.47(17) and 0.63(23) MeV 
respectively; other energies changed by less than 1%. For A = 8 nuclei, the GFMC improves 
the '^T energy by ~10 MeV; our best '^v provides only ~1.5 MeV of this improvement. 

Table |V1II| shows excitation energies for the A = 8 nuclei. The errors are the combined 
statistical errors of the energies computed for the given state and the ground state. We see 
that, while the absolute VMC energies have substantial errors as compared to the GFMC 
values, the VMC excitation energies are typically within an MeV of the GFMC values. 



Figure ^ shows the energies of nuclear states for 4 < A < 8 and Fig. jTU] shows the 
corresponding excitation spectra. As is the case for the ground-state values, most of the 
A = Q and 7 excited-state energies have been recomputed with improved VMC wave func- 
tions and all of the GFMC propagations have been made to r = 0.2 MeV~^. The VMC 
excitation energies of the lowest ^Li(i ) and ^Li(| ) states have been significantly reduced 
from the values reported in Ref. These are the result of improved mixings of the different 
symmetry states in the variational wave function. 

The figures confirm the principal A = Q and 7 results of Ref. for A = 8: the Hamil- 
tonian consisting of Argonne f is plus Urbana IX underbinds nuclei in the p-shell with the 
underbinding becoming worse as one increases A or N — Z. However the predictions of the 
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excitation spectra are generally reasonable; the rms error of all of the 16 excitation energies 
computed by GFMC is only 650 keV; for the 9 states with experimental width less than 
200 keV it is 540 keV. Considering that no parameters were adjusted to fit these p-shell 
energies, this is quite respectable. We do note that, where there is an experimental value to 
compare to, the VMC excitation energies of second excited states of a given J'^, which are 
very difficult to compute by GFMC, are usually too high. 

Although the Argonne f is plus Urbana IX Hamiltonian consistently underbinds the p- 
shell nuclei, the errors are small compared to the magnitudes of the potential energies. Table 
IX| shows the perturbatively extrapolated GFMC expectation values for the kinetic energy 
and potential energy terms. As is discussed in Ref. the perturbatively extrapolated 
GFMC terms do not add up to the total energy, which is the most reliable number computed 
by this method. The total two-body potential, Vij is dominated by the one-pion exchange 
term, vfj. There is a large cancellation between the kinetic and two-body potential energies so 
that their sum is only 15% to 20% of the two-body potential. The total three-body potential 
is typically 4% to 5% of the two-body potential. The repulsive V^^^ typically cancels 45% 
of the V^jl- As an indication of the smallness of the errors produced by this Hamiltonian, 
consider the ground state of ^Be: the difference of the computed and experimental energies 
is 2.1 MeV; the two-body potential is -301 MeV, V^l = -27 MeV, and V^f^ = +12.3 MeV. 
Thus the corrections that have to be made to Vijk are significantly smaller than the terms 
already included. 



B. Isobaric Analog States 

Energy differences of isobaric analog states are sensitive probes of the charge- 
independence-breaking parts of the Hamiltonian. To study these it is useful to express 
the energies in an isobaric multiplet, characterized by A and T, in terms of the isospin 
multipole operators of order n: 

Ea,t{T.) = E «S!TQn(T,T,) . (5.2) 

n<2T 

The Qn{T,Tz) are orthogonal functions for projecting out isovector, isotensor, and higher- 
order terms the first terms are Qo = 1, Qi = Tz, and Q2 = |[3T| — T(T + 1)]. The 
coefficients a^"^ are then obtained from the calculated energies: 
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E Q„(T, Tz)Ea,t{T.)/ E Ql^T, T,) , (5.3) 



or perturbatively from expectation values of the isomultipole operators present in the Hamil- 
tonian. 

In first-order perturbation theory, the electromagnetic interaction contributes to the a*^") 
for n = 0, 1, and 2, the kinetic energy to n = and 1, the nuclear CSB potential to 
n = 1, and the nuclear CD potential to n = 2. Because a significant portion of the vf^ 
comes from one-pion exchange, there should also be a CD component to V^^l^ however, a 
plausible extension of the Urbana IX model gives negligible contributions of 3 keV or less 
to the n = 2 terms. The a*^") for higher n are zero in first order with our Hamiltonian, and 



there is little experimental evidence for n > 3 terms in nuclei |3^. We have made VMC 
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calculations of the a^^^ and a^"^^ in first order by using a CI wave function of good isospin, T, 
and simply varying to compute the Ea,t{Tz):, in the GFMC calculations we evaluate 
expectations of the isomultiplet operators for the = — T nuclei. The CSB and CD parts 
of the Hamiltonian can induce corresponding changes in the nuclear wave functions, leading 
to higher-order perturbative corrections to the splitting of isospin mass multiplets. However, 
it is difficult for us to estimate these higher order effects reliably in either the VMC or GFMC 
calculations. 

Table pCl shows results for the T = 1 and 2 isovector and isotensor coefficients in A = 8 

— (2) 

nuclei compared to experiment. In obtaining the experimental a\ [ coefficient, we have used 
the average of the two isospin- mixed (2"*"; 0+1) states; our calculation of this mixing is dis- 
cussed in the next subsection. The contributions of the complete f v^^ , v'"^^, and K'~'^^ 
terms to the a*^"-' are also given. The results show that the present Hamiltonian and VMC 
wave function match the experimental CSB and CD of the T = 1 multiplet fairly well, but 
the calculated value for 03^2 is about 10% too small. The GFMC consistently reduces the 
coefficients, mostly by reducing the Coulomb term. For the isovector coefficients, this wors- 
ens the comparison with experiment, however the case of a^^\ is significantly improved. In 
all these calculations, the Coulomb interaction between protons is the dominant contribu- 
tion to the a*^"'-', but the strong interaction CSB terms serve to bring the final results closer 
to experiment. However, in view of the fact that the A = 8 nuclei are underbound with 
the present Hamiltonian, it is premature to use these calculations as a precise test of the 
charge-independence-breaking components of the interaction. 



C. Isospin Mixed States 

The (2"'";1) isobaric analog of the ^Li ground state in '^Be is very close in energy to the 
second (2"'";0) excitation, as shown in Fig. 0, and is in fact experimentally observed to be 
isospin- mixed with it. There are also fairly close (1"^;0,1) and (3"*"; 0,1) pairs at slightly higher 
energies in the *Be spectrum. Our VMC calculations do not get the (2"*";0,1) states so close 
together, although the other pairs do come out quite near each other. However, we can 
calculate the isospin-mixing matrix elements that connect these pairs of states: 

Eoi(J) = (vI/(J+,0)|i/|vI/(J+,l)) . (5.4) 

This is done by increasing the model space to include both T = and 1 components, 
which corresponds to a wave function vector with 10,752 terms, and employing the same 
off-diagonal evaluations used to determine the l3Ls[n] wave function components. 

Results for the Eqi{J) are given in Table |X|. The experimental values are determined 



from the observed decay widths and energies ||3^. The dominant contribution, from the 
Coulomb potential, typically accounts for less than half of the matrix element. We find the 
magnetic moment part of the electromagnetic interaction and the strong CSB interaction 
can provide a significant boost, although we still underpredict the experimental mixing by 
~20%. The spatial symmetry components of the different wave functions, as given in the 
Tables pT| and for the (2+;0,l) and (3+;0,l) pairs are indeed fairly close to each other, 
since they have similar sizes and signs for the largest components. However, the (1"'";0,1) 
states are not so similar, particularly due to the large T = 1 ^P[31] component that is not 
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available to the T = state. This may be why the -£"01(1) is noticeably smaller, which shows 
up through the change in sign of the magnetic moment contribution. 



VI. MOMENTS AND DENSITY DISTRIBUTIONS 

The proton rms radii, magnetic moments, and quadrupole moments for the A = 8 
nuclear ground states are given in Table |X11|. The calculations have been made in impulse 



approximation for both the initial VMC wave function and perturbatively after GFMC 
constrained path propagation. Comparing the VMC and GFMC radii, it appears that the 
variational wave functions are a little too compact for this Hamiltonian. However, given 
that the present interaction underbinds these nuclei, the actual radii should be smaller 
than the GFMC values. Attempts have been made to determine the matter radius of ®He 
by proton-scattering experiments in inverse kinematics. Interpretation of the data is not 
model- independent, however, and has resulted in estimates ranging from 2.45 fm to 2.6 
fm The corresponding VMC and GFMC matter radii from Table PCII| are 2.69 fm and 



2.92 fm, respectively. In the not too distant future, it may be possible to trap ^He, *Li, and 
long enough to determine their charge radii to high precision by atomic means. 
The magnetic moments of ^Li and have been measured by /3-radiation detection 
of implanted polarized ions [^. Our lA calculations should be supplemented by meson- 
exchange currents (MEG); their contributions are dominantly isovector, and have been shown 
to change the magnetic moments of '^H and ^He by ± 0.4 a? |^ . The isoscalar average of the 



calculated magnetic moments is close the experimental value, and it is plausible that MEG 
contributions will bring both magnetic moments into good agreement with experiment. The 
quadrupole moments of ^Li and ^B have been measured by /?-NMR and nuclear quadrupole 
resonance techniques ||38|. The VMC and GFMC results are fairly close for ®Li, and just a 
little above the experimental value. The calculations for ^B are significantly further apart, 
despite the proton rms radii differing by only a small amount; nevertheless they bracket the 
experimental value. We have also calculated the hexadecapole moment for these two ground 
states, and find it to be consistent with 0. 

One- and two-nucleon density distributions have been calculated in both VMC and 
GFMC. The GFMC densities are a little more spread out and less peaked than the cor- 
responding VMC densities, as reflected in the charge radii differences; here we show only 
the GFMC densities. Single-nucleon density distributions for the A = 8 nuclei are shown 
in Fig. |TT]; they are normalized such that the integrated value equals the appropriate total 
value of or Z. The two protons in '^He are the most peaked distribution, which should be 
expected on the grounds that they are mostly confined to the a core, while the six neutrons 
in ^He have the broadest distribution. The proton distribution in ^Li is also fairly peaked 
near the origin, but is broader than in ®He since there is one additional proton in the p-shell. 
The neutron distribution in ^Li is comparable to the ^He neutrons near the origin, is slightly 
larger in the 1-2 fm range, and then falls below at larger distances; the intermediate-range 
excess may be due to the significantly greater binding of ^Li compared to ^He. In contrast, 
the ^Be proton and neutron distributions are much less peaked at the origin and are rather 
fiat out past 1 fm. This could be because ^Be has a significant 2a component in its intrinsic 
structure, with the two a's sitting side-by-side; as discussed in the next section. 

A logarithmic plot of the single-nucleon densities for ^'^'^He isotopes is shown in Fig. O. 
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The neutron halos in ^He and ^He are clearly evident, while the proton cores of these nuclei 
are nearly identical. The peak neutron and proton distributions in ^He and ^He are much 
reduced compared to "^He, because of the motion of the a core against the center-of-mass. 
Two-proton density distributions for ^'^'^He are shown in Fig. |13[ The relative proton-proton 
density in ^'^He is not effected by the motion of the a relative to the center of mass; thus ppp 
is an indicator of changes in the internal structure of the a core in ^'^He. Here we note that 
there is only a small change in the a core of ^'^He, i.e., about a 10% reduction in the peak 
compared to ^He, which must be compensated in the long-range tail, since the total integral 
under the curves is unity in each case. This reduction could be due to swelling of the a core, 
or due to charge-exchange interactions between the protons and p-shell neutrons. 

We have also calculated the single-nucleon momentum distributions shown in Fig. 
using the VMC wave functions and the method of Ref. |^9|], with some algorithmic im- 
provements. The N{k) are normalized to N/{2'kY or Z/(27r)'^. All the distributions show a 
remarkably similar structure: a low-momentum core attributable to s- and p-shell orbitals, 
followed by a high-momentum tail beyond 2 fm~^ that is the sum of many small-amplitude 
higher-orbital contributions. This high-momentum tail is already evident in the ^He mo- 
mentum distribution and is intimately related to the D-state of ^He, present due to the 
strong tensor forces and corresponding correlations in the wave function. The ^He proton 
distribution is very similar to that of "^He, showing that the a core is not much altered, while 
the neutrons exhibit a p-shell peak near 0.5 fm~^. The ^Li proton distribution is a little 
broader, as might be expected by the addition of a p-shell proton, while the neutrons show 
a smaller p-shell peak. The ®Be momentum distribution is almost exactly double that of 
■^He except at very low momenta, indicating a significant 2a structure. 



VII. INTRINSIC SHAPES 

The ground state and first two excited states of '^Be have an approximate rotational 
energy spectrum, and are assumed to be well approximated as two a's rotating around their 
common center of mass. This structure is not manifest in the shell model part of the VMC 
wave functions, in which the ^Be states are constructed from an a-like core surrounded by 
four p-shell nucleons coupled to the appropriate total angular momentum. In this section 
we describe an attempt to study the intrinsic structure of the ^Be states described by the 
VMC wave functions including correlations. 

The standard Monte Carlo method for computing one-body densities, p(r), is to make a 
random walk that samples |\Ef (ri, r2, ■ ■ ■ , Ta) P and to bin ri, r2, ■ ■ ■ , for each configuration 
in the walk. The density is then proportional to the number of samples in each bin. In the 
case of a J = nucleus, this "laboratory" density will necessarily be spherically symmetric. 

We can attempt to find the intrinsic density in body-fixed coordinates by computing the 
moment of inertia matrix: 



A 

m=y: 

i=i 



ViXi yf ViZi I , (7.1) 



for each configuration. We then find the eigenvalues and eigenvectors of M., rotate to 
those principal axes, and bin the resulting r'^, r'g, ■ ■ ■ , r^. The eigenvector with the largest 
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eigenvalue is chosen as the z' axis; further choices of the ± direction along z' and which 
eigenvectors to assign to y' and x' may also be made or averaged over. This procedure will 
not produce a spherically symmetric distribution, even if there is no underlying deformed 
structure, because almost every random configuration will have principal axes of different 
lengths and the rotation will always orient the longest principal axis in the z' direction. 
We have made a number of tests using simple wave functions with no internal correlations 
and also nuclei like ^He which, in our models, should have no intrinsic deformations. We 
find that the projected "intrinsic" density for such cases is always prolate but that no other 
artificial structure is introduced. 

When the above procedure is applied to the '^Be rotational states, a dramatic intrinsic 
structure is revealed. Figure |15| shows this calculation for the ground state. The figure shows 
contours of constant density plotted in cylindrical coordinates, with the z-axis being the axis 
of quantization. The calculation was made for the VMC wave function. The left side of 
the figure shows the standard density calculation; we can think of this as the density in the 
"laboratory" frame. For the J = ground state, this is spherically symmetric as shown. 
The right side of the figure shows the intrinsic density computed as described above. In this 
case the orientation (up or down) along z' and around z' (r = + y^) was averaged over. 
It is clear that the intrinsic density has two peaks, with the neck between them having only 
one-third the peak density; we regard these as two a's. 

Figure |l^ shows the corresponding calculation for the J = 4+, Mj = 4, state of ^Be. In 
this case the the laboratory density does not have to be spherically symmetric and there 
is evidence of two a's rotating around the z axis in the equatorial plane. The projection 
to the intrinsic frame rotates these to the z' axis and results in an intrinsic density that 
is insignificantly different from that obtained for the ground state. A calculation for the 
J = 2+ state also produces the same intrinsic density. 

These results, obtained for the VMC wave functions, suggest that the 0^, 2"'", and 4"'" 
wave functions for ^Be have the structure of a deformed rotor consisting of two a's. The 
structure is not manifest in the shell model part of the VMC wave functions; it is induced 
by the correlations. The optimum spatial correlations between pairs of s-shell and p-shell 
nucleons are similar, but those between one s-shell and one p-shell nucleon are different. 
Calculations using GFMC configurations give very similar results except that, especially for 
the 7 = 4"*" state, the rms radius is still growing at the end of the GFMC propagation. 

If the O"*", 2+, and 4"*" states are generated by rotations of a common deformed structure, 
then their electromagnetic moments and transition strengths should all be related to the 



intrinsic moments. Table pCIIl| shows the VMC computed values of the quadrupole and 
hexadecapole moments, and of the B(E2) and B(E4) strengths for these states. In the 
rotational model these are related by simple Clebsch Gordon coefficients to the intrinsic 
quadrupole, Qo, or M4 moment depending on the multipolarity, A. These moments are 
defined as 



Qo = y dr'p{r'y'Y2o{r') , (7.2) 
and 

M4 = I rfr'p(r')r'V4o(r') , (7.3) 
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where p is the point proton density and r' refers to the intrinsic (body-fixed) frame. The 
last column gives the extracted values of these intrinsic moments; we see that, with the 
exception M4 for the J = 2+ state, these are remarkably constant. 

We can also compute values of Qq and M4 by integrating over the projected body-fixed 
densities. These values might be too large because, as described above, the projection 
method can introduce an excessive deformation in the intrinsic shape. The values for the 
0+, 2+, and 4+ states are Qo = 26.2, 27.9, and 26.7, and M4 = 55, 62, 64; the Qo are in good 
agreement with the values given in the Table, which were obtained from the spectroscopic 
quadrupole moments. We note that the ratios of these Qq and M4 values are in reasonable 
agreement with the ratio obtained for a diatomic model of ^Be assuming two point a's 
separated by 4 fm. 

We can attempt to project out an intrinsic structure of the T = 0, J = 1+ and J = 3^, 
states of *Be. In this case we also find two peaks but of somewhat smaller density than for 
the 0+, 2+, and 4+ states. Clearly these peaks cannot be due to two a's rotating around 
each other. They presumably refiect occurrence oi a + t + p and a-l-^He-l-n structures in 
these states. 

VIII. CONCLUSIONS 

We have made quantum Monte Carlo calculations for the ground states and low-lying 
excitations oi A = 8 nuclei interacting by realistic two- and three-nucleon potentials. These 
calculations have been made practical by the development of a constrained-path algorithm 
for the complex spatial and spin-isospin wave functions needed to describe these nuclei. This 
algorithm greatly reduces the "fermion-sign problem" that quantum Monte Carlo methods 
are subject to, allowing us to obtain binding energies for a given Hamiltonian that are 
accurate to 1 to 2%. 

The Hamiltonian we have used, consisting of the Argonne vis two-nucleon and Urbana 
IX three-nucleon potentials, gets the general features of the light p-shell nuclei fairly well, 

inchiding the bulk of the experimentally observed binding and the correct ordering and 
approximate spacing of the excitation spectrum. However, with accurate calculations, we 
can identify specific failings of this Hamiltonian, including a few percent underbinding in 
N — Z nuclei, which gets progressively worse as the neutron-proton asymmetry increases, 
and spin-orbit splittings in the excitation spectrum that are too small. Nevertheless, we can 
see that the energy differences with experiment are much smaller than the magnitude of the 
short-range (and least-well-known) part of our three-nucleon interaction, so it is plausible 
that new three-nucleon force models may be able to reproduce the ground and excited states 
near the 1% level. 

We have also been able to study energy differences between isobaric analog states and 
the isospin-mixing matrix elements in ^Be. While final conclusions should be reserved until 
the bulk energies are corrected by improved three-nucleon potentials, we can see that the 
charge- dependent and charge-symmetry-breaking components of the Argonne Vig potential 
are making significant contributions that improve the agreement with experiment. 

Finally, we have also studied the moments, densities, momentum distributions, and in- 
trinsic shapes of these nuclei. It appears that the observed static moments may be under- 
stood with our present microscopic model after we include meson-exchange current contri- 
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butions. We also see evidence for strong clustering in the light p-shell nuclei, particularly 
the 2a character in ^Be. This clustering is built into the variational wave functions by the 
strong pair correlations, which depend upon the nucleon orbits and seem to be preserved in 
the GFMC propagation. 

With the present method, and the continuing rapid increases in computational power 
that massively parallel machines are bringing, we are confident that we will be able to extend 
our calculations to the ^4 = 9, 10 nuclei in the near future, and to ^^C in a few year's time. 
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TABLES 



TABLE I. Values of variational parameters in p-shell nuclei; in the expression for Vp, L is the 



total orbital angular momentum. 










1'3L[31] 


i'=^L[22] 




3l[211] 








0.85 




0.85 




0.85 




0.85 


dsp 


(fm) 




3.2 




3.2 




3.2 




3.2 


Cpp 






0.1 




0.2 




0.3 




0.4 


dpp 


(fm) 




3.2 




3.2 




3.2 




3.2 


Vp (MeV) 




-20.0+L 




-20.0+L 




-20.0+L 




-20.0+L 


Rp (fm) 




4.0 




4.0 




4.0 




4.0 


Qp (fm) 




1.5 




1.5 




1.5 




1.5 




TABLE IL PlsIu] components for ^ 


He states, listed in order of increasing 


I excitation. 


(J-T) 




1 


S[22] 




iD[22] 






3p[211] 


(0+;2) 






0.842 










-0.539 


(2+;2) 










0.958 






-0.287 


(i+;2) 
















1 


(0+;2) 






0.546 










0.838 


(2+;2) 










0.294 






0.956 




TABLE in. 


PLS[n] components for 


^Li states, listed in order of increasinj 


5 excitation. 


(J^ 


T) 


=^P[31] 


=^D[31] 


=^F[31] 


ip[31] 


iD[31] 


iF[31] 


3S[22] 


3D[22] 


(2+ 


,1) 


0.936 


-0.289 


-0.121 




0.118 






-0.104 


(1+ 


1) 


0.740 


-0.334 




0.576 






-0.067 


-0.068 


(0+ 


1) 


1 
















(3+ 


1) 




0.912 


0.275 






-0.296 




-0.067 


(2+ 


,1) 


0.343 


0.854 


0.386 




0.046 






0.047 


(1+ 


1) 


0.644 


0.410 




-0.618 






-0.162 


-0.091 


(1+ 


1) 


-0.036 


0.812 




0.531 






-0.229 


0.063 


(4+ 


1) 






1 












(2+ 


1) 


-0.111 


-0.260 


0.550 




0.780 






-0.101 


(3+ 


,1) 




-0.269 


0.922 






0.037 




-0.275 
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l/VJDijIli i V . PIjS\ 


„] components for ^Be states, listed in order of increasing excitation. 




T) iS[.l] 




■^P[31] 




■■^F[31] 


(0+ 


0) 0.997 




-0.076 






(2+ 


0) 


0.999 


0.033 


0.037 


-0.018 


(4+ 


0) 


0.997 






0.079 


(2+ 


0) 


-0.019 


0.948 


-0.314 


0.045 


(1+ 


0) 




0.861 


0.508 




(1+ 


0) 




-0.526 


0.850 




(3+ 


0) 






0.891 


0.453 


(4+ 


0) 


-0.089 






0.996 


(2+ 


0) 


-0.038 


0.270 


0.844 


0.463 


(3+ 


0) 






-0.457 


0.889 


(0+ 


0) 0.098 




0.995 








TABLE V. ^Li(gs) results (in MeV) using 


different tlircc-body propagators. 




Propagator 


/Y2n,A) 


/y2^,Cv 






G 




O.Ul 1 






-31.25(12) 


G,n 


= 1 


-8.5(2) 


-5.2(1) 




-31.15(11) 


G,n 


= 4 


-8.9(2) 


-5.3(1) 




-31.23(7) 


G,n 


= oo 


-8.8(2) 


-5.2(1) 




-31.25(12) 


just two-body 


-5.8(2) 


-3.9(1) 




-29.86(23) 






liVBLE VI. Scaling (^f calcula 


ti(ni with tiVHlcuu size 








A 


P = A{A-\)I2 


2^ X 7(A, T) 




n(x^Be) 


^He 


4 


6 


16x2 




0.001 


^Hc 


5 


10 


32x5 




0.010 


6Li 


6 


15 


64x5 




0.036 


^Li 


7 


21 


128x14 




0.33 


8Be 


8 


28 


256x14 




1 


8He 


8 


28 


256x20 




1.4 


8Li 


8 


28 


256x28 




2 


9Be 


9 


36 


512x42 




8.7 




10 


45 


1024x42 




24 


iiB 


11 


55 


2048x132 




200 


12c 


12 


66 


4096x132 




530 


8n 


8 


28 


256x1 




0.071 


14^ 


14 


91 


16384x1 




26 
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TABLE VII. Experimental and quantum Monte Carlo energies of A=2-8 nuclear ground states 



in MeV. 


Ay{ TIT. rp\ 


V iViv^ y^J^ J 


V iviL^ y j 




tliX.pi 




00/18/' 






00/1 R 
— Z.ZZ4D 


3-fT(l + . 1\ 


— 0.10(^1 j 


— O.OZ(^i j 


Q A'7/-\\ 


Q /IQ 
— 0.4o 




— ZD.y7(oj 


— z7.7o(oj 


— zo.34(4j 


00 on 
— ZO.30 


^HefO+- 1) 

X A \ \J J J- J 


-23 64f7l 


-24.87(7) 


-28 llf9) 


-29.27 


6Li(l+;0) 


-27.10(7) 


-27.83(5) 


-31.15(11) 


-31.99 


^He(i-;i) 


-18.05(11) 


-19.75(12) 


-25.79(16) 


-28.82 


^Li(i";^) 


-31.92(11) 


-33.04(7) 


-37.78(14) 


-39.24 


«He(0+;2) 


-17.98(8) 


-19.31(12) 


-27.16(16) 


-31.41 


8Li(2+;l) 


-28.00(14) 


-29.76(13) 


-38.01(19) 


-41.28 


8Be(0+;0) 


-45.47(16) 


-46.79(19) 


-54.44(19) 


-56.50 
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TABLE VIII. VMC, GFMC, and experimental excitation energies (adjusted to their respective 
ground states) in MeV. Numbers in parentliescs are Monte Carlo error estimates for theory, and 
uncertainties in the energy for experiment. We also give experimental widths where known in keV 
(except where otherwise noted). The ^Be states marked with an * are isospin-mixed. 





T) 


VMC 


GFMC 


Expt 


r 






(McV) 


(MoV) 


(McV) 


(keV) 


8Hc(2+ 


:2) 


2.22(16) 


3.01(23) 


3.59 




8He(l+ 


;2) 


3.32(17) 


4.42(25) 






«He(0+ 


;2) 


4.72(18) 








«He(2+ 


;2) 


5.01(17) 








8Li(l+; 


1) 


1.34(18) 


0.57(24) 


0.98 


12(4) fs 


8Li(0+; 


1) 


2.83(19) 


1.91(29) 






8Li(3+; 


1) 


3.35(18) 


2.66(28) 


2.26 


33(6) 


8Li(2+ 


1) 


3.86(18) 








8Li(l+, 


1) 


4.22(19) 




3.21 


K. 1000 


8Li(l+; 


1) 


5.32(19) 




5.4 


K 650 


8Li(4+; 


1) 


6.08(18) 


6.27(31) 


6.53 


35(15) 


8Li(2+; 


1) 


6.20(18) 








8Li(3+: 


1) 


7.31(18) 




6.1 


« 1000 


8Li(0+; 


2) 


11.24(18) 


12.10(25) 


10.82 


< 12 


8Be(2+ 


0) 


2.39(25) 


2.91(25) 


3.04(3) 


1500(20) 


«Be(4+ 


0) 


9.95(24) 


9.58(27) 


11.4(3) 


3500 


8Be(2+ 


1) 


18.60(23) 


18.02(27) 


16.63* 


108(1) 


8Bc(2+ 


0) 


20.29(24) 




16.92* 


74(1) 


8Be(l+ 


1) 


19.89(23) 




17.64 


11(1) 


8Be(l+ 


0) 


20.03(23) 


18.09(33) 


18.15 


138(6) 


8Be(l+ 


0) 


21.73(24) 








8Bc(3+ 


1) 


21.77(22) 




19.07 


270(20) 


8Be(3+ 


0) 


21.85(24) 


19.53(33) 


19.24 


230(30) 


8Be(4+ 


0) 


25.85(22) 




19.86 


700(100) 


8Be(2+ 


0) 


23.61(23) 




20.1 


^ 1100 


8Be(3+ 


0) 


25.53(23) 








8Be(0+ 


0) 


28.62(26) 




20.2 


< 1000 


8Be(0+ 


2) 


29.55(22) 


29.93(25) 


27.49 


5.5(2.0) 
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TABLE IX. Kinetic and potential energy contributions to GFMC energies in MeV. 





K 








13 






2H(1+;0) 


19.81 


-22.04 




0.018 


-21.28 








51.2(3) 


-58.9(3) 


-1.2(0) 


0.04(0) 


-45.0(2) 


-2.2(0) 


1.0(0) 


%e(0+;0) 


110.7(7) 


-135.3(7) 


-6.3(1) 


0.86(0) 


-102.4(5) 


-11.5(1) 


5.2(1) 


6He(0+;l) 


138.(1) 


-164.(1) 


-7.0(1) 


0.86(0) 


-119.(1) 


-13.1(2) 


6.1(2) 


6Li(l+;0) 


153.(1) 


-183.(1) 


-7.3(2) 


1.68(1) 


-143.(1) 


-13.6(2) 


6.3(2) 


^He(i-;i) 


156.(2) 


-183.(2) 


-8.1(2) 


0.87(0) 


-131.(1) 


-15.2(4) 


7.2(3) 




190.(2) 


-227.(2) 


-9.1(2) 


1.77(1) 


-172.(1) 


-17.9(3) 


8.8(2) 


8He(0+;2) 


166.(1) 


-195.(1) 


-8.0(2) 


0.87(0) 


-135.(1) 


-15.8(2) 


7.8(2) 


8Li(2+;l) 


217.(2) 


-257.(2) 


-10.2(2) 


1.86(1) 


-192.(1) 


-20.8(4) 


10.6(3) 


«Be(0+;0) 


248.(2) 


-301.(2) 


-14.9(3) 


3.27(1) 


-227.(1) 


-27.1(4) 


12.3(3) 



TABLE X. VMC and GFMC isovector and isotensor energy coefficients a~2rr ^'S^- 



(n) 




j^CSB 






yCSB 


+ ^CD 




Total 




VMC GFMC 


VMC 


GFMC 


VMC 


GFMC 


VMC 


GFMC 


Expt. 


"8,1 


24 


20 


1689 


1586 


66 


63 


1779(4) 


1669(7) 


1770 


"8,1 








155 


132 


-11 


3 


144(4) 


135(7) 


145 


"8,2 


17 


15 


1478 


1423 


40 


47 


1535(3) 


1485(5) 


1659 


"8,2 








140 


129 


26 


29 


166(4) 


158(4) 


153 






TABLE XL 


VMC 


isospin-mixin 


g matrix elements for ^Be in 


keV 








yCoul 


^MM 


yCSB 




^CSB 




Expt. 


2+ 




63 


18 


26 




2 


109(4) 


149 


1+ 




40 


-1 


15 




1 


55(2) 


120 


3+ 




33 


15 


13 




1 


62(2) 


63 



TABLE XIL VMC and GFMC values for point proton rms radii (in fm) , for magnetic moments 
(in /xjv), and quadrupole moments (in fm^) in impulse approximation. In the GFMC calculations, 



the rms radii of the two excited states of 


Be are 


steadily 


growing 


with propaj 


Ration time. 






1/2 










Q 






VMC 


GFMC 


VMC 


GFMC 


Expt. 


VMC 


GFMC 


Expt. VMC 


8Hc(0+;2) 


1.91(1) 


1.98(1) 














8Li(2+;l) 


2.21(2) 


2.18(2) 


1.12(1) 


1.3(1) 


1.653 


3.2(1) 


3.6(2) 


3.19(7) 0.1(2) 


«Be(0+;0) 


2.25(1) 


2.48(1) 














8b(2+;1) 


2.47(3) 


2.57(1) 


1.46(1) 


1.4(2) 


1.036 


5.5(2) 


7.9(4) 


6.83(21) 0.7(8) 


8C(0+;2) 


2.90(4) 


3.17(1) 
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TABLE XIII. VMC calculation of moments and transition strengths for different multipolari- 
ties, A, in the rotational states of ^Be; the units are fm and the proton charge. 



Observable 


A 


Value 


Qo or Mi 


Q(2+) 


2 


-7.61(8) 


26.6(3) 


Q(4+) 


2 


-9.76(10) 


26.8(3) 


B(E2, 2+ ^ 0+) 


2 


14.8(4) 


27.3(4) 


B(E2, 4+ ^ 2+) 


2 


18.2(4) 


25.3(3) 


(2+,2|r4y4,o|2+,2) 


4 


2.5(2) 


53.(4) 


(4+,4|r4y4,o|4+,4) 


4 


4.3(2) 


34.(2) 


B(E4, 4+ ^ 0+) 


4 


167.(16) 


39.(2) 
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FIGURES 

FIG. 1. The experimental spectrum for ^ = 8 nuclei. All narrow low-lying states are shown, 
as well as the wide rotational states in ^Be. Higher-lying and negative parity states in ^Be and 
some wider states in ^Li are not shown. 

FIG. 2. Constrained and unconstrained propagation tests for ^Li(gs). See the text for a de- 
scription of the symbols. 

FIG. 3. Constrained and unconstrained propagation tests for an eight-neutron drop. See the 
text for a description of the symbols. 

FIG. 4. Constrained and unconstrained propagation tests for ^He(gs). See the text for a de- 
scription of the symbols. 

FIG. 5. Constrained propagation to large r for states of ^Li. 

FIG. 6. GFMC E{t) for states of ^Be. The solid lines show the averages of the last 7 values; 
the dashed lines show the corresponding statistical errors. 



FIG. 7. GFMC E{t) for states of ^Li. 

FIG. 8. GFMC E{t) for states of ^He. 

FIG. 9. VMC, GFMC, and experimental energies of nuclear states for 4 < A < 8. The light 
shading shows the Monte Carlo statistical errors or experimental uncertainties. 

FIG. 10. VMC, GFMC, and experimental excitation energies of nuclear states for 6 < ^ < 8 

FIG. 11. The neutron and proton densities in ^He, ^Li, and ^Be. 

FIG. 12. The neutron and proton densities in ^He, ^He, and ^He, shown on a logarithmic scale. 

FIG. 13. The proton-proton densities in ^He, ^He, and ^He. 

FIG. 14. The neutron and proton momentum distributions in ^He, *Li, and ^Be. 
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FIG. 15. Contours of constant density, plotted in cylindrical coordinates, for ^36(0"*"). The left 
side is in the "laboratory" frame while the right side is in the intrinsic frame. 

FIG. 16. Contours of constant density for ^Be(4+). 
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Fig. 1 (Wiringa, et al.) 
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Fig. 2 (Wiringa, et al.) 
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